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Abstract 

In this paper, the algorithm for determining the 
stencil of a one- dimensional Essentially Nonoscillatory 
(ENO) reconstruction scheme on a uniform grid is 
reinterpreted as being based on extrapolation. This 
view leads to another extension of ENO reconstruc- 
tion schemes to two-dimensional unstructured triangu- 
lar meshes. The key idea here is to select several cells 
of the stencil in one step based on extrapolation, rather 
than one cell at a time. Numerical experiments confirm 
that the new scheme yields sharp nonoscillatory recon- 
structions and that it is about five times faster than 
previous schemes. 

Introduction 

ENO (Essentially Nonoscillatory) schemes devel- 
oped by Harten, Osher, Shu and others [1,2, 3, 4] are a 
class of higher order numerical schemes for solving hy- 
perbolic partial differential equations. These equations 
frequently admit solutions in which some smooth struc- 
ture is interspersed with discontinuities. ENO schemes 
are designed so as to achieve a high order of accuracy in 
the smooth regions and at the same time avoid spurious 
oscillations at discontinuities. In such cases, conven- 
tional fixed stencil schemes either damp out the smooth 
structure while obtaining a monotone shock profile or 
capture the smooth structure and incur spurious oscil- 
lations near the shock. 

The key feature of ENO schemes is the use of an 
adaptive stencil in the reconstruction step, i.e. in con- 
structing point values of u(x) from the set of cell 
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averages fij. For each cell, a searching algorithm de- 
termines which stencil of neighboring points has the 
smoothest data and u(z) is then reconstructed using 
data only from this stencil. The variable stencil has 
the effect of avoiding 0(1) spurious oscillations near 
discontinuities. 

With unstructured grids offering greater geometric 
flexibility and ease of adaptation, many efforts [5,6, 7, 8] 
have been made to extend ENO schemes to unstruc- 
tured grids. Harten and Chakravarthy [5] describe a 
cell centered extension of the finite volume ENO scheme 
to triangulated unstructured grids along with details 
of implementation although no numerical results are 
presented. * Abgrall [6,7] presents an ENO cell vertex 
scheme where the control volumes are the cells of the 
dual mesh and presents results for ID convection and 
gas dynamics. 

In this paper, we first reinterpret the basic one di- 
mensional ENO scheme on a uniform grid as a scheme 
based upon extrapolation. Based on this point of view, 
we present another extension of ENO schemes to un- 
structured meshes. This new scheme has the advantage 
that it is more efficient than the scheme of Harten and 
Chakravarthy[5] and gives almost identical results on 
fine enough meshes. We consider only reconstruction 
in this paper. Implementation of this scheme in a flow 
solver will be the subject of future work. 

ID ENO Reconstruction Revisited 

Consider m - 1-th order ENO reconstruction of 
a function u(x) using the primitive function technique 
(see [1] for details). Let the cell of interest be [xj t Xj+i] 
with cell average tij+ 1 / 2 * The reconstruction will use 
the cell averages [iij+i/ 2 -m+ii —«j+i/ 2 +m-i]* The 
primitive function of u(x) is defined by 

h(x)= Turn (i) 

Jx 0 
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and is known at the cell faces up to an additive con- 
stant. An m-th order interpolant Q(x) that interpo- 
lates H{x) at [x,-, ...x t + m ] is constructed next and the 
reconstruction, denoted by Vj+\| 2 {x)^ is then obtained 
as 

f;+i/2(*) = ^ (2) 

There are precisely m choices for i namely [j — m -f 
1, ... .j]. The nonoscillatory properties of the recon- 
struction result from choosing i so that [x,-, ...Xj+ m ] is 
the smoothest possible region around the cell [xj y x ; -+i]. 
One way of doing this is by a recursive procedure, as 
described below. 

We begin with i = j, the linear interpolant for 
Q(x) for which the stencil is [xj, Xj+i], and keep adding 
cells until we have the required number. Assume that 
at some stage in this process the stencil is [xj, ....x,+ p ]. 
At the next step, the choice is between the two stencils 
[x„ ...JCi+P+d or [X|_ i , . ..x t + p ]. If H[xi t .. .Xi+p+i] is less 
than i/[x;_i, ...Xj +p ], where H[..] denotes the divided 
differences of H{x), the former stencil is picked, oth- 
erwise the latter. In current practice, it is customary 
to weight the two divided differences being compared 
so as to favor stencil movement towards the centered 
stencil. However, this is not crucial to the equivalence 
we wish to establish below and thus has been dropped 
for clarity. 

Let us now consider a different method of sten- 
cil selection. We begin as before and assume that at 
some stage in this process the stencil is [x, , ...Xj+ P ]. We 
reconstruct tt(x) based on this stencil and denote this 
reconstruction by u p (x). Now, we extrapolate v p (x) to 
the neighboring cell [xj+ p , x»+ p +i], evaluate the cell av- 
erage of i^(x) over this cell, and subtract from this the 
true cell average of this cell. Let us denote this resid- 
ual as R+ . The same thing can be done for the other 
neighboring cell [x*_i,Xj] and we end up with another 
residual which we denote as iZ”. If |iZ + | < |iZ”|, we 
choose [xi } ...x,*+p+i] for the stencil and proceed to the 
next step. If not, the stencil chosen is [x,_i, ...x,*+ p ]. 

Since the reconstruction at any stage is the deriva- 
tive of a primitive function, the values of the residuals 
can be calculated as 

p - 1 

= H[xit £*+ P +i] n (*£»+?+ i ^t'+o (3) 

fc =0 

P-1 

R- = H [ar, x { +p ] (*<+*+, - ar<_i) (4) 

k = 0 

It follows that for a uniform grid this method of sten- 
cil selection is identical to the standard ENO method 


via divided differences of H since the two products in 
Eqs. (3) and (4) are the same. As in standard ENO, 
the residuals can be appropriately weighted so as to 
favor stencil movement towards the centered stencil. 
In smooth regions, the residuals are 0(Ax p ) while if 
the extrapolated cell contains a jump discontinuity the 
residual is 0(1), entirely analogous to scaled divided 
differences. For nonuniform grids in one dimension, 
the two methods differ because of different weighting 
factors on the divided differences and it is not clear to 
us which is the better choice. Limited numerical exper- 
iments with the new method yield solutions that are 
similar to that from standard ENO. 

The extrapolation method outlined above is rem- 
iniscent of subcell resolution introduced by Harten[9]. 
There, the extrapolation is performed to locate the dis- 
continuity within a cell. Unlike sub cell resolution, the 
extrapolation method can be readily extended to multi- 
dimensions as we show below. In addition, ENO inter- 
polation can also be viewed as being based on extrap- 
olation. In this case, the extrapolated point value is 
compared to the true point value, their difference be- 
ing the appropriate residual. 

In what follows below, we apply this extrapolation 
idea to the problem of ENO reconstruction on a two- 
dimensional unstructured mesh. 

ENO Reconstruction on Unstructured Grids. 

An unstructured triangular mesh is characterised 
by a set of vertices (x^y*), where i = 1, ..nv and a set 
of connection coefficients c(i y j ) where i = 1, nc and j — 
1,3. nc is the number of cells or triangles and nv the 
number of vertices. For the i-th. triangle, c(i,j) specify 
the labels of its vertices. A given set of vertices may be 
connected in many ways to form an unstructured grid. 
A systematic way of connecting a given set of points is 
the so called Delaunay triangulation which has several 
desirable properties. 

We will consider only cell centered schemes in this 
work so that the control volumes or cells will be the tri- 
angles themselves. The reconstruction problem can be 
stated as follows: Given the cell averages of a piecewise 
smooth function tx(x), reconstruct the function within 
a particular cell in a nonoscillatory manner to within a 
specified order of accuracy. For some theoretical results 
on this problem, the reader is referred to Abgrall[6,7]. 

The first step is to consider the case where u(x) is 
globally smooth. Let 0» be the cell of interest whose 
centroid is chosen as the origin of the coordinates (x, y). 
For an m-th order reconstruction we seek the recon- 
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structed function as 

m 

v(x,y) = '%2 Y2 a h kx3 y k ( 5 ) 

1=0 j+k=l 

Since there are (m + l)(m + 2)/2 unknowns, it is clear 
that we need a stencil of at least that many cells near 
Ci. For conservation, one of these cells must be C*. 
How these cells are selected lies at the heart of the ENO 
reconstruction procedure and will be discussed below. 
For now, assume a stencil of cells with cell indices Jq{i), 
where q = 1, ...(m4- l)(m + 2)/2 have been found. The 
linear system for determining can be written as 

Yh Y2 a i- k < x3yk > i = ( 6 ) 

1=0 j + k=l 

where q = 1, ...(m 4* l)(m + 2)/2 and 

<xV > ? =T 7 ^-r // (7) 

|W ff l JJc Jq 

There is no guarantee that the linear system (6) will 
be nonsingular, but when it is, the solution of it yields 
the required reconstruction. It is clear, however, that 
the matrix will be singular only for special cases. 

We now describe a method of selecting the stencil 
due to Harten and Chakravarthy[5] which we shall refer 
to as ENOHC. We begin with the single triangle Ci for 
which the zero-th order reconstruction is ao,o = or 

v(x) = Ui (8) 

Let us denote by Jp{i) a stencil of p cells at the p-th 
stage of the stencil selection process where 

J p (i) = **p} ( 9 ) 

These p triangles can be used to fit p terms of the expan- 
sion^) so that associated with this stencil are p values 
of While these p terms can be taken in any order, 
it is natural to take them in the lexicographic order, 
namely, 1, x, y, x 2 , xy } y 2 , .... Let J * denote cells not in 
J p (i) that share a common side with cells in J p (i). In 
the next step we add to J p {i) any one triangle from J* 
and fit the p 4- 1 terms in the expansion (5) using the 
p + 1 triangles. We do this for all triangles in J* and 
pick the one that minimizes the sum of the absolute val- 
ues of the p+ 1 coefficients a. The stencil now contains 
p + 1 triangles and p 4- 1 values of and we proceed 
to the next step. When p = (m 4- l)(m 4- 2)/2 — 1 
the process is complete and we have the required m-th 
order ENO reconstruction. 

The cost of adding one cell to the stencil can be 
roughly estimated as follows. For every element of J*, 


calculating the norm of the coefficients a involves cal- 
culating 2 p 4 ’ 1 quadratures (the additional row and 
column of the linear system ( 6 )) one p 4-1 x p + 1 ma ~ 
trix inversion and a norm calculation. Neglecting the 
cost of finding the minimal element, the approximate 
cost of adding one cell to the stencil may be written as 

|j;|((2p+l)Q + M(p+l) + iV(p+l)) (10) 

where Q is the cost of a quadrature, M(n) is the cost 
of inverting an n x n matrix and N(n) is the cost of a 
norm calculation of n terms. 

This method seems to be a logical extension of the 
standard ENO scheme in one dimension. Numerical 
tests that we have performed with the method (to be 
presented below) show that sharp nonoscillatory recon- 
structions are obtained away from the boundary. A 
curious feature of the scheme is that the stencil chosen 
and hence the reconstruction obtained are not coordi- 
nate invariant, which is due to the fact that terms are 
fitted one by one. The main drawback of the scheme 
is its expense, especially at higher orders. This may 
be traced to the fact that the number of triangles re- 
quired is proportional to m 2 and triangles are added 
one by one. To go from an m-th order reconstruction 
to m 4 - 1 -th order requires the addition of m 4 - 2 cells to 
the stencil, each of which is added in a separate step. 
In the next section, we introduce a simplification that 
will allow us to choose many cells together resulting in 
a more efficient algorithm. 

Stencil Selection via Extrapolation 

As before, we begin with the single triangle C, 
and the zero-th order reconstruction given by Eq.( 8 ). 
Assume we have an l-th order reconstruction for the 
cell of interest. This means we have a stencil J p (i ) of 
p - (/ 4 - 1 )(/ 4 - 2)/2 cells and a reconstruction of the 
form 

v <( x ) = Y^ S a j, kX *y k ( n ) 

q = 0 j+k = q 

As before, we have the set of neighbors to this stencil 
J * . For every candidate cell in q G we compute the 
residual given by 

= f v,{x,y)dxdy ( 12 ) 

\^q\ J JC q 

If the cell Cq and all cells in the stencil J p (i) lie in 

a smooth region, the residual is 0(h l + l ) where h is 

some representative cell size. If not, the residual is 
0(1). Hence, we compute this residual for all candidate 
triangles and pick those / 4 - 2 triangles that have the 
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least value of the residual. Once we select the / + 2 
additional triangles, we perform a linear reconstruction 
problem to obtain the / + 1-th order reconstruction. 
Then we proceed to the next stage until when / = m 
we have the desired reconstruction. 

The key difference between ENOHC and this 
method is that we select all / + 2 triangles in one step. 
This is also the main problem with this method which 
works best when every cell has a large number of neigh- 
bors. In many cases, especially for the lower orders, 
there are not enough triangles. For example, while go- 
ing from the zero-th order reconstruction to a first order 
reconstruction, we have to select 2 triangles from the 3 
neighbors. This is usually not enough of a choice and 
the reconstruction is often oscillatory. Numerical tests 
suggest two further modifications which are necessary 
for a robust reconstruction algorithm. The first is we 
have found it necessary to limit the order of the recon- 
struction based on the magnitude of aj Specifically, 
we compute the norm 


m = 


i 

(/ + l)(/ + 2)/2 


E E k*i 


g-0 j+k=q 


(13) 


and if N(l) > 1/ft we terminate the reconstruction at 
that stage. Unlike ENOHC, this has the effect of re- 
ducing the order of the reconstruction near discontinu- 
ities, in the same spirit as Suresh[10]. There does not 
seem to be a clear choice for the mesh refinement pa- 
rameter h for unstructured meshes. In our numerical 
experiments, h was chosen to be the minimum of the 
circumcircle radii. 

The second modification is based on the detectabil- 
ity of the residuals. As mentioned above, the residual 
for a cell in the smooth region is 0(h* +1 ) while for a cell 
containing a discontinuity it is 0(1). Thus detection of 
the “bad” cells is most difficult at the lower values of 
/. On the other hand, ENOHC is relatively inexpensive 
at the lower orders. Hence, the first few stages can be 
obtained using ENOHC and the latter stages obtained 
using extrapolation. Thus the second modification that 
we have used is to obtain a first order reconstruction us- 
ing ENOHC and use the extrapolation procedure from 
then on. In numerical tests below, we shall refer to the 
scheme described above as ENOEX. 

We remark that a scheme in between ENOHC and 
ENOEX can also be constructed which like ENOHC 
adds cells one by one but uses extrapolation to select 
the cells. Here, at the p-th stage of stencil selection, we 
compute the residual (12) for all cells in J * and pick a 
single cell that has the least residual. This triangle is 
added to J p and the associated linear system solved to 


get the new p - hi values of aj t * and we proceed to the 
next stage. The effort involved in going from a stencil of 
p cells to p+ 1 cells is one quadrature for every element 
of J* and one matrix inversion or 

\j;\Q + M(p+l). (14) 

This is still a substantial saving over Eq.(10) although 
much more expensive than ENOEX. We do not consider 
this intermediate scheme any further. 

Numerical Experiments 


We describe below several numerical tests per- 
formed on the new ENOEX reconstruction algorithm. 
The unstructured grids used for these purposes were 
generated by applying random perturbations to the ver- 
tices of a structured grid and then applying the Bowyer- 
Watson Delaunay triangulation algorithm to the result. 
Typical grids are shown in Fig. 1 and Fig. 9, where 
the domain of reconstruction is [—1, 1] x [—1, 1]. For the 
coarse grid (Fig. 1) a structured grid of 40 x 40 was used 
and the magnitude of the random perturbation was set 
to ±40% of the structured mesh size (i.e. 0.05). For the 
fine grid, a structured grid of 80 x 80 was used with the 
same relative magnitude of the perturbations. The per- 
turbations are set to zero on the boundary. The cell size 
parameter in ENOEX was chosen to be the minimum 
of the circumcircle radii over the whole grid. This is 
obtained from the Bowyer- Watson algorithm and fixed 
apriori in the computations presented below. 

In studying the accuracy of the reconstructions, 
contour plots of functions on unstructured meshes can 
be made with presently available software tools. How- 
ever, getting one-dimensional line plots of the function 
along curves is often difficult and involves further in- 
terpolation. For this reason, we have found it conve- 
nient to perform the reconstruction on to the nodes of a 
suitably chosen structured mesh. From the structured 
mesh values it is quite easy to study one-dimensional 
line plots along the coordinate directions. The struc- 
tured mesh is chosen so that there is roughly one node 
per cell and is 41 x 41 for the coarse grid and 81 x 81 
for the fine mesh. 

1 . Reconstruction of a smooth function 

The first tests are performed on smooth functions 
and serve as a convenient debugging mechanism. We 
reconstructed the smooth function 

u(x,y) = sin(7r(x + y)) 

on the two unstructured meshes shown in Figs. 1 and 
9 respectively using ENOEX4 (fourth order ENOEX 
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reconstruction). The errors in the max-norm and com- 
puted order of accuracy are shown in Table- 1. These 
errors were computed at the centroid of each cell. Since 
there is no well defined mesh refinement parameter, we 
have calculated the order of accuracy based on four dif- 
ferent parameters, namely the maximum and minimum 
circumcircle radii and the square root of the maximum 
and minimum cell areas. The numerical order of accu- 
racy varies from 6.28 to 3.55, depending on the mesh 
refinement parameter used. In theory, all should ap- 
proach 5 on fine enough meshes. 

2. Reconstruction of a piecewise smooth function 

Since we deal only with the reconstruction prob- 
lem, the initial function to be reconstructed has to be 
sufficiently general. A suitable function has been con- 
structed by Abgrall [6] which has several curved discon- 
tinuities and smooth structures. The function (which 
is slightly different from that in Abgrall [6]) is defined 
as follows: If <j) is any angle and r = — (y — xtan(<} i>)/3 
then f<j>(x,y) is defined as: 


r if 
i if 
if 


r < -1/3, 4(x, y) = — r sin(?rr 2 /2), " 

r > 1/3, 4(x, y) = 2r - sin(37rr)/6, > 

M < 1/3, f*(x,y) = |sin(27rr)| , 


(15) 


and the function to be reconstructed u(x } y) is 


u(x, y) = y) ^ x < cos(7ry/2) 

«(*, y) = y) + cos(2iry) (16) 

if x > cos(7ry/2) 


case, the two reconstructions are almost identical with 
the only difference being at about (—0.4, —0.9) where 
ENOEX4 is a bit noisy. 

The CPU time comparisons are presented in Table 
2. Since a lot of time in these runs is spent in com- 
puting the intersections with the structured grid and 
in the computation of the cell averages, these times 
were not included in the run times listed in Table-2. It 
can be seen that ENOEX4 scheme is about five times 
faster than ENOHC4. These comparisons are meant to 
serve only as a guide as they depend on the efficiency of 
programming. Furthermore, the gap narrows for lower 
orders being identical for m = 1 . 

Conclusions and Remarks 


In this paper, we have proposed a new essentially 
nonoscillatory (ENO) reconstruction procedure on two- 
dimensional unstructured meshes. At the heart of the 
method is a reinterpretation of standard ENO recon- 
structions as being based on extrapolation. Numerical 
experiments confirm that the reconstruction is indeed 
nonoscillatory and that the new fourth order scheme is 
about five times faster than earlier schemes. 

We view ENOEX as a first step towards faster 
ENO schemes on unstructured meshes. Many issues 
need to be investigated further. Weighting of the cen- 
tral stencil, which is important for stability consider- 
ations, can be accomplished in several ways, with and 
without some notion of distance from the cell of inter- 
est. Better choices of the cell size parameter may exist. 
These questions will be explored further in future work. 


The cell averages of this function over each cell are com- 
puted by using an adaptive quadrature routine(from 
the IMSL library), where in most cells, the estimate of 
absolute error is machine zero. However, for cells near 
discontinuities, the absolute errors are larger, with the 
largest error being O(10 -6 ). 

Figures 3 and 4 show the results of reconstruction 
using ENOHC4 and ENOEX4 (fourth order reconstruc- 
tions) respectively on the coarse mesh shown in Figure 
1. Note that this reconstruction corresponds to a fifth 
order scheme for convection or gas dynamics. Figure 
2 is a contour plot of the exact solution. It can be 
seen that the results of ENOEX4 are comparable to 
ENOHC4 although the latter has slightly better reso- 
lution at discontinuities. Figures 5-8 show line plots 
of the reconstructions at various x and y sections of 
the domain. For the most part, the reconstructions are 
similar with ENOEX showing a slightly larger under- 
shoot in Figure 7. Figures 11 and 12 show correspond- 
ing results for the fine mesh shown in Figure 9. In this 
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Table - 1 


Grid 

Max Error 

Cmax 


V A m in 

m 

Coarse 

9.8698 xlO -4 

0.0542 

0.0155 

0.0100 

0.0566 

Fine 

3.9099 xlO" 5 

2.8077(-2) 

6.2788(-3) 

6.0046(-3) 

0.0304 

Order 

— 

4.90 

3.55 

6.2870 

5.1962 


Table 1: Reconstruction of the smooth function sin(7r(x + y)) on un- 
structured meshes. C min and C max are the minimum and maximum of the 
circumcircle radii and A mtn and A max are the minimum and maximum cell 
areas respectively. 


Table - 2 


CPU Times (Ymp Seconds) | 

Grid 

ENOEX4 

ENOHC4 

Coarse 

21.694 

124.252 

Fine 

100.07 

525.121 


Table 2: Comparison of CPU times for fourth order reconstruction using 
ENOHC and ENOEX. 
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Figure 1: Unstructured Coarse Grid. 1600 vertices and 3024 
triangles. Generated by random perturbations of a 40 X 40 
uniform structured grid. 


Figure 2: Exact solution on the coarse grid 
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Figure 9: Unstructured Fine Grid. 6400 vertices and 12482 p 1Q; Exact so , ution Qn the fine grid 

triangles. Generated by random perturbations of an 80 X 
80 uniform structured grid. 
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a * EN0EX4 
O = EN0HC4 
— = Exoct Solution 
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